function [V,T] = delaunay_lai(x,y)
% This function returns a Chebyshev Triangulation of the unit square.
% This matlab program is copyrighted @2001 by Ming-Jun Lai and Paul Wenston
% through University of Georgia Research Foundation, Inc..
% Find a bounding box
xmin=min(x); xmax=max(x);
ymin=min(y); ymax=max(y);
xdiff=xmax-xmin; ydiff=ymax-ymin;
x1=xmin-xdiff; x2=xmax+xdiff;
y1=ymin-ydiff; y2=ymax+ydiff;
p=[x1 x2 x2 x1; y1 y1 y2 y2];
t=[1 1; 2 3; 3 4; 0 0];
[m,n]=size(x);
if m> n
x=x'; y=y';
end;
% Triangulate
[p,t]=delaun(p,t,x,y);
% Remove the  bounding box
p(:,1:4)=[];
d=find((t(1,:)>0&t(1,:)<5)|(t(2,:)>0&t(2,:)<5)|(t(3,:)>0&t(3,:)<5));
t(:,d)=[];
t(4,:)=[];
t(1:3,:)=t(1:3,:)-4;
T=[p(:,t(1,:));p(:,t(2,:));p(:,t(3,:))];
T=T';
T = add_element(T,x',y');
V = [x',y'];

function [p,t,c]=delaun(p,t,c,x,y)
%  [p,t,c]=delaun(p,t,c,x,y) 
%This matlab program constructs Delaunay triangulation from a list of vertices
%whose x-components in x and y-component in y with bounding box p and inital
%triangles in t. 
if nargin==4
  y=x;
  x=c;
  c=circumc(p,t);
end;
np=size(p,2);
tp=size(t,2);
p=[p,[x;y]];
t=[t,zeros(4,2*length(x))];
c=[c,zeros(3,2*length(x))];
scale=max(max(abs(p)));
small=100*eps;
tol=small*scale;
for l=1:length(x)
  i=find((c(1,1:tp)-x(l)).^2+(c(2,1:tp)-y(l)).^2-c(3,1:tp).^2<=-tol);
  M=sparse(t([1 2 3],i),t([2 3 1],i),ones(3,length(i)),np,np);
  M=M-M'>0;
  [j,k]=find(M);
  i2=length(j);
  tp=tp+2;
  i1=[i tp-1 tp];
  t(1,i1)=j';
  t(2,i1)=k';
  t(3,i1)=ones(1,i2)*np+1;
  c(:,i1)=circumc(p,t(:,i1));
  np=np+1;
end;

function c=circumc(p,t)
%This matlab program computes the circumcircle of a triangle with vertices and
%triangles 
x1=p(1,t(1,:)); x2=p(1,t(2,:)); x3=p(1,t(3,:));
y1=p(2,t(1,:)); y2=p(2,t(2,:)); y3=p(2,t(3,:));
a11=x2-x1; a12=y2-y1; b1=x2.^2+y2.^2-x1.^2-y1.^2;
a21=x3-x1; a22=y3-y1; b2=x3.^2+y3.^2-x1.^2-y1.^2;
d=2*(a11.*a22-a21.*a12);
cx=(b1.*a22-b2.*a12)./d;
cy=(a11.*b2-a21.*b1)./d;
c=[cx;cy;sqrt((x1-cx).^2+(y1-cy).^2)];

function NewT = add_element(T,X,Y)
  n = size(T,1);
  NewT = [];
  for i = 1:n
    v1 = find(X==T(i,1) & Y==T(i,2));
    v2 = find(X==T(i,3) & Y==T(i,4));
    v3 = find(X==T(i,5) & Y==T(i,6));
    NewT = [NewT;v1,v2,v3];
  end